#include <iostream>
#include <vector>

#include <map>
#include <set>

#include <cassert>

#include <sds/sfile.h>
#include <sds/sdataset.h>
#include <sds/sdataview.h>

#include <matrix.h>
#include <cluster.h>
#include <random.h>

using namespace std;
using namespace sds;

void renumber_dfs(int i, int j, int n, int c, 
		 SDataView<unsigned char>& col,
		 SDataView<unsigned int>& num) {
  if (i < 0 || j < 0 || 
      i >= col.getSize()[0] || j >= col.getSize()[1] ||
      col[i][j] != c) return;
  num[i][j] = n;
  col[i][j] = 255;
  renumber_dfs(i - 1, j, n, c, col, num);
  renumber_dfs(i + 1, j, n, c, col, num);
  renumber_dfs(i, j - 1, n, c, col, num);
  renumber_dfs(i, j + 1, n, c, col, num);
}

void renumber(SDataView<unsigned char>& col,
	      SDataView<unsigned int>& num, int& n) {
  n = 1;
  for (int i = 0; i < col.getSize()[0]; ++i) {
    for (int j = 0; j < col.getSize()[1]; ++j) {
      if (col[i][j] == 255) continue;
      if (col[i][j] != 0) {
	renumber_dfs(i, j, n, col[i][j], col, num);
	++n;
      } else {
	num[i][j] = 0;
      }
    }
  }
}

void Read(const string& filename, const std::string& segname,
	  SDataView<Point>& src, SDataView<unsigned int>& seg,
	  vector<SInteger>& size, int& segcnt) {

  SFile segfile(segname);
  SDataView<unsigned char> segcol = 
    segfile.getDataset(0).getDataView<unsigned char>();
  size = segcol.getSize();
  
  seg = SDataView<unsigned int>(size);

  renumber(segcol, seg, segcnt);

  SDataView<unsigned char> source[CH];
  SFile file(filename);
  for (int k = 0; k < CH; ++k) {
    SDataset dataset = file.getDataset(k);
    source[k] = dataset.getDataView<unsigned char>();
  }

  src = SDataView<Point>(size);
  for (int i = 0; i < size[0]; ++i) {
    for (int j = 0; j < size[1]; ++j) {
      Point point;
      for (int k = 0; k < CH; ++k) {
	point[k] = source[k][i][j];
      }
      src[i][j] = point;
    }
  }  
}


int main(int argc, const char *argv[]) {
  assert(argc == 3);

  SDataView<unsigned int> seg;
  SDataView<Point> src;
  vector<SInteger> size;

  int n;

  Read(argv[1], argv[2], src, seg, size, n);

  --n;

  vector<Cluster> segments(n);

  for (int i = 0; i < size[0]; ++i) {
    for (int j = 0; j < size[1]; ++j) {
      if (seg[i][j] != 0) {
	segments[seg[i][j] - 1] += src[i][j];
      }
    }
  }

  for (int k = 0; k < n; ++k) {
    Matrix covariance(CH, CH);
    for (int i = 0; i < CH; ++i) {
      for (int j = 0; j < CH; ++j) {
	covariance(i, j) = segments[k].covariance(i, j);
      }
    }
    Vector eval = covariance.eigenValues();
    for (int i = 0; i < CH; ++i) {
      std::cout << eval(i) << '\t';
    }
    std::cout << std::endl;
  }


  return 0;
}
